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This paper presents ANSYS Fluent simulation results and analysis for self-pressurization 
of a flightweight, cryogenic, liquid hydrogen tank in 1-g. These results are compared with 
experimental data, in particular, pressure evolution and temperature measurements at a set 
of sensors. The simulations can be analyzed to identify and quantify heat flows in the tank. 
Heat flows change over time and influence the self-pressurization process. The initial rate of 
self-pressurization is sensitive to the initial temperature profile near the interface. 
Uncertainty in saturation pressure data and the accuracy of experimental measurements 
complicate simulation of self-pressurization. Numerical issues encountered, and _ their 
resolution, are also explained. 


Nomenclature 
A = cross-sectional area, perpendicular to heat flux, of tank wall or tank fluid level, m 
C,C,C, = Modified Lockheed equation constants 
CC, = Specific heat (solid), specific heat at constant pressure (fluid), J/kg-K 
C, Cli Cvap = NOF fraction, liquid, vapor at a point, unitless, [0, 1] 


D, = degradation factor 

k = thermal conductivity, W/m-K 

g = acceleration due to gravity, m/s” 

GH2, LH2= hydrogen gas, liquid hydrogen 

L = length scale of VOF interface thickness, m 

MLI = Multi-Layer Insulation 

MW = molecular weight, kg/kmol 

Metux: m = mass flux due to evaporation (>0) and condensation (<0), kg/s-m’, volumetric rate, kg/s-m°* 
f = MLI layer density, cm 

N = number of MLI layers 

NIST = National Institute of Standards and Technology 

D, Psa(T) = static pressure, Pa, saturation pressure at temperature, Pa 

P = MLL interstitial gas pressure, 1.x10° torr 

Q = heat flow, W 

dup = heat flux through wall, through MLI blanket, W/m? 

shroud = net radiation heat flux to shroud, W/m? 

Ry = universal gas constant, J/kmol-K 

Teulk = Temperature of the bulk of liquid, K 

T, T,,(P) = static temperature, saturation temperature at pressure, K 

Ty, TTayg = MLI blanket temperatures, hot edge, cold edge, and their average, K 

Tinterface = Temperature of the liquid/vapor interface, K 

UDF = User Defined Function (Fluent) 

VOF = Volume of Fluid 
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Z = Compressibility factor, real gas molar volume to ideal gas molar volume, unitless 
a = underrelaxation parameter, unitless 

é = radiative emissivity, unitless 

viscosity, Pa s 

p density, kg/m* 

Ocond Sevap= Mass accommodation coefficients for condensation and evaporation, unitless 

lig, vap, boil = subscripts indicating liquid and vapor phases, and boiling conditions 


I. Introduction 


ASA is interested in long-term, in-space storage of 

cryogenic propellants to support future exploration 
missions, including upper stages and potentially propellant Chania ii 
depots. Cryogenic propellants promise higher specific impulse —_"*n#"" 4 
than storable hypergolic fuels, but storage for long duration 
missions must be demonstrated. Improving the capabilities of 
computational tools to predict fluid dynamic and 
thermodynamic behavior in cryogenic propellant tanks under 
settled and unsettled conditions is research supported under 
NASA’s Evolvable Cryogenics project. 

To test these simulation tools, the focus of this work is self- 
pressurization experiments performed in the early 1990’s at 
NASA Lewis research center with a flightweight liquid 
hydrogen tank. The tank was tested at the heat fluxes, and 
conditions expected in future space exploration missions. The 
tank had high performance multilayer thermal insulation and a Figure 1: Cryogenic hydrogen tank without 
low mass-to-volume ratio. The intent of these original tests MUI insulation. From NASA TN D-8320. 
was to understand the underlying physics of 
self-pressurization and predict design 
performance. 

This paper is organized by topic and issue 
; =a (sub-section) to cover the many issues in this 
Z x EY complex simulation. 


II. Experimental Geometry and Grid 


The cryogenic hydrogen tank simulated is 
from the K-Site experiment (Figure 1) 
performed at the Cryogenic Propellant Tank 
Facility at NASA Plumbrook research center 
in the 1990’s. This large, flightweight, 2219 
aluminum tank with MLI insulation was tested 
in vacuum conditions [1] at 350 K_ shroud 
temperature and 1-g. The 2.2 m diameter tank 
consists of two elliptical domes and a very 
small barrel section. Although the tank was 
tested for self-pressurization at 29%, 49%, and 
83% fill levels, the current simulations are 
focused on the 49% fill level. MLI 
performance is documented in [2]. 

The computational domain is _ two- 


Wall |.) 


Witt +1 | A : . : . . : 
i dimensional, axisymmetric. The baseline grid 


a + 
Figure 2: Cryogenic hydrogen tank grid (left), with inset 
(top) of unstructured grid for the lid, and inset (right) 
showing the grid where the tank wall meets the liquid/vapor 
interface. A thicker wall band with transition to thinner, 
chemically milled walls is also shown. The lid is treated as solid. 
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contains 14,800 cells. A multi-block 
structured grid (Figure 2) captures variations in 
wall thickness: chemically milled wall sections 
(2 mm thick) (Figure 1), thicker wall bands 
(4.2 mm thick) where support struts are 


attached. Further, an unstructured grid represents the complex geometry of the lid of the manhole access. Some 
grids have a vent boundary condition to solve for the initial condition. 

Solution accuracy is not simply proportional to the numerical scheme’s order of accuracy, but also the local 
solution gradient; grid stretching has a role. Compared to the bulk fluid in the tank interior, this solution gradient is 
1000 times greater near the walls (momentum and thermal boundary layers) and the interface (thermal layer). 
Consequently, the baseline grid has 0.5 mm normal spacing at the walls and 0.5 mm spacing at the liquid/vapor 
interface—S0 times more than the coarsest grid spacing in the bulk fluid region. 

Coarser, similar, and finer grids were used to assess grid convergence (Table 3); however the finer grid is 
prohibitively expensive for long, time-dependent simulations. 


III. Material Properties and Numerical Methods 


ANSYS Fluent version 16.0 [3] is used to solve thermal equations in the solid walls coupled to thermal/fluid 
equations in the fluid region (two-dimensional, axisymmetric, Navier-Stokes equations). The simulation includes the 
VOF equations to capture two-phase flow, and the fluid flow is modeled as laminar flow. This transient simulation 
is second-order in space and time. 


A. Material Properties 

Constant fluid physical properties at reference conditions, (20.354 K, 103.632 kPa), do not adequately capture 
thermal conductivity, k, in the gas phase and in the aluminum tank walls over the simulation’s temperature range. 
Linear equations in temperature accurately represent NIST data [4] for viscosity, 4, and thermal conductivity, k, of 
the gas phase hydrogen. NIST data [5] for 5083 aluminum is used for piecewise linear representations of the 
specific heat, C, and thermal conductivity, k. 5083 aluminum data is indistinguishable from 2219 aluminum data [6] 
for thermal conductivity, and NIST data is publicly available and readily accessible on the internet [4] [5]. 
a real 
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Hydrogen, 
gas, differs from an 
ideal gas by up to 12% 
over the temperatures 
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se “18-20 and pressures in the 

a ~ 16-18 simulation, as shown 

= “1416 in Figure 3. The 
12-14 


compressibility factor, 
Z, deviates most near 
the saturation line, at 


“10-12 
=8-10 


Saturation Line 


Temperature 
(kK) 


2 a high pressures, and at 
215 654 low temperatures. In 
21 aa. the simulations, the 
20.5 ideal gas assumption is 
20 used, so Z is not 


corrected. 

It is important to 
Figure 3: Percentage deviation of a real gas from ideal gas from NIST [4] note that the ratio of 
hydrogen data: (1-Z)*100%. Positive values indicate reduced volume. The black thermal conductivities 
curve is the saturation line, P,,,(T). Note the temperature scale change at 25 K. for GH2:LH2:Al is 

1:6:1180—the heat 
flux through the metal walls is much, much higher than in the vapor or liquid, and only partially mitigated by thin 
walls. Also, note that the volumetric heat capacity, oC,, is as much as 40 times higher for LH2 than GH2, hence the 
liquid can absorb much more heat than vapor for the same temperature rise. Aluminum’s volumetric heat capacity, 
pC, is similar to that of GH2. 
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B.Saturation Conditions 
Why is the saturation line, P,a(T), and its representation so important? The self-pressurization experiment’s 
vapor phase pressure evolution is the focus of this work. Vapor pressure is essentially constant in the gas phase at 
any time—perturbations travel everywhere at the speed of sound. Further, vapor pressure is not only the saturation 
pressure, Pea(Tinterface), at the liquid/vapor interface, but it reflects the temperature of the interface, Tinterface- 
Deviations from a single value of Tinerjace—cold or hot, liquid or gas, impinging on the interface—results in 
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evaporation in warm areas and condensation in cold areas that pushes the interface temperature to a single interface 
value—all mediated by the vapor pressure being held virtually constant at any time throughout the vapor phase. 
Hence, accurate saturation conditions are important if pressure evolution is to be captured. 

Accurate representation of the saturation line is also important because saturation pressure is very sensitive to 


temperature—in fact, 
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Figure 4: Saturation line, P,q(T), by NIST [4] and 
Reynolds [8] data for parahydrogen. Visually, the cubic 
interpolant lies on and obscures the NIST and Reynolds data. 
Uncertainty bars (40.2%) are plotted for NIST data, but they 
are not visible on this scale. The Clausius-Clapeyron 
equation is not accurate over the temperature range. 
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three multiplications and three additions. Experimentally [1], the 
initial condition is Pyq(20.354 K) = 103,632 kPa. 

Fluent mentions the Clausius-Clapeyron equation [9, p. 608], 
and although it is accurate very close to one temperature, it is 
inappropriate for representing saturation pressure over large 
temperature ranges (Figure 4). 

Clearly, the experimental temperature accuracy makes it 
difficult to match simulation results with experimental self- 
pressurization data. 


C. Mass Transfer: Evaporation and Condensation at the 
Interface 

The rate of mass transfer is calculated from the Schrage 
equation, Eq. (1), which is derived, in turn, from Maxwell’s 
distribution. This evaporation/condensation internal boundary 
condition enforces saturation conditions at the liquid/vapor 
interface. Eq. (1) is implemented in a Fluent User-Defined- 


~30 kPa/K over the range of conditions studied here. 


Further, the accuracy of 


experimental liquid-vapor temperature 
measurements (silicon diode transducers) is +0.1 
K [1, p. 2], and wall temperatures are accurate to 
+0.6 K. Consequently, the accuracy of 
temperature measurements lead to +3 kPa 
accuracy of pressure measurement. Ref. [7] gives 
a good explanation of silicon diode temperature 
sensors, their calibration and operation. 

The saturation line conditions for hydrogen, 
Psa(T), are represented by a cubic polynomial 
least squares fit to NIST data [4]; the interface 
temperature range is [20.25, 23.5] K, 
corresponding to gas phase pressures of 97.188 to 
230.123 kPa, as shown in Figure 4. This 
polynomial fits NIST data to within +7 Pa across 
this temperature range which is less than the 
estimated uncertainty of 0.2% (~200 Pa). With 
Horner’s rule, the cubic can be implemented with 
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Figure 5: Experimentally observed effect of 

starting conditions on pressure rise rate. An 

~8 kPa pressure difference develops in the first 

two hours. g = 3.5 W/m’, Fill Level 83%, 


From [1] and [22]. 


Function (UDF) with the approximation Tijg = Tyap, Sevap = Ocona and calculated every time step—but not every 


iteration. 
m = 2 MWyap (« Psat(T liq) = @ Pyap ) (1) 
f lux 2G 2nRy evap Tig cond Tay > 
; m : 
net = ~~ where L=1/|Vc| converts flux to volumetric rate (2) 
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The length scale of the VOF interface thickness, L, converts the flux to a mass transfer rate [10], Eq. (2). To 
smooth this value in time, it is underrelaxed every time step, 1, Mfreigx = amr. + (1 — a) mez}, where the 
underrelaxation factor, a, is typically 0.02. 


I. Mass Accommodation Coefficients: 

Although there is evidence that mass accommodation coefficients are near 1, here Oeyayp = Ocong = 0.001 is used. 
An accommodation coefficient near 1 corresponds to a very fast mass transfer rate. We believe that another rate- 
limiting step exists in the process. In particular, heat flow through thin thermal boundary layers at the liquid/vapor 
interface limits the heat available for evaporation/condensation. Some mass transfer schemes are based on the 
interfacial heat flux [11]. The current simulations are not sensitive to the accommodation coefficient value, o. 

Classical molecular dynamics computer simulations of the air/water interface give mass and thermal 
accommodation coefficients of 0.99 and 1.0, respectively, at 300K [12]. Experimentally, the thermal 
accommodation coefficient was shown to be 0.84+0.05 in argon at 271K using sound resonance in a spherical 
chamber [13]. The experiment creates strong sound waves with a pressure oscillation both above and below 
saturation pressure—so evaporation and condensation alternate with the pressure wave; a liquid layer can exist on 
the walls but without a net heat flow [14]. 


2. Unphysical Conditions and Loss of Fluid Mass: 

Condensation in an all liquid cell—or evaporation in an all vapor cell—must be avoided since it leads to a loss of 
fluid mass (see Section VI-A) when using Fluent’s Define_Mass_Transfer() UDF. At every iteration and for every 
cell, Eq. (3) is applied. Multiplication by cjjg or Cygy ensures zero evaporation or condensation, respectively, in 
unphysical conditions, and this is Fluent’s suggested approach [9, p. 608]. 


if (Mreiax = 0) {mh = 0} 
else if (Myeiax > 0) {rn = Mretax * Ctig } Evaporation—liquid must be present (3) 
else if (Mreiax < 0) {rm = Mretax * Cia Condensation—vapor must be present 


Adding Eq. (3) to the mass transfer scheme avoids mass loss, but it may make convergence more difficult. For 
simulations that are otherwise identical, slightly shorter time steps have been required in simulations using Eq. (3). 


D. Heat Transfer Boundary Condition on the Tank Exterior Surface 
Heat transfer on the tank wall exterior surface is a boundary condition that includes heat leaks through MLI 
insulation and through struts, plumbing and 


oe hs wires plus radiation to/from the shroud. 

= SS The insulation consists of two MLI blankets 

Re ~/\V\S® each with 15 layers of double aluminized 

oS Mylar film (Figure 6) plus reinforced Mylar 

= face sheets, with silk net spacers between 

BS these 17 layers. Table 1 gives the 

= experimentally observed heat leaks where 

: aS total heat flow is 48.2 W. The heat leaks 

LH2/GH2 Tank Inner MLI Outer MLI = Radiation shroyg through struts, plumbing and wires are 
Wall Blanket Blanket in Vacuum located on the wall, scaled to Table 1 


values, and added to the MLI heat flux. 


Figure 6: Architecture of MLI insulation layers. See [2, pp. 8-9]. The local heat flux, guyz, through the 


Not to scale. 


Table 1: Experimentally measured heat fluxes for Tprouq = 350 K, P’ = 4.0e-6 torr [21, p. 8], 


Heat Leak Contributions Heat Leak (W) % Total Heat Leak Heat Flux (W/m’) 
12 Struts 2.813 5.83 0.201 
Plumbing & Ducts 3.194 6.62 0.228 
Wires 0.879 1.82 0.063 
MLI 41.352 85.72 2.952 
Total Heat Leaks 48.239 100. 3.444 
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MLI insulation layers surrounding the tank is modeled with the Modified Lockheed equation [15], Eq. (4), and 
radiation exchange with the shroud, qshroua, Eq. (5). These heat fluxes are calculated each time step for each cell and 
applied on the tank wall exterior surface as a thermal boundary condition. A Fluent UDF uses Eq. (4) and the local 
wall temperature to find the local heat flux, gjy;;, through each MLI blanket—but it must be the same through each 
blanket and equal to the external radiative heat flux, Gshroug at each grid cell on the external wall. The algorithm [16] 
that enforces Gyr inner =QMLI outer = Yshroua finds MLI blanket edge temperatures, 7, and T,, for equal heat flux 
through each blanket based on the local tank wall temperature, plus external radiation based on shroud temperature. 


dur = Dg[C(0.017 + 7 x 10-°(800 — Tayg) + 2.28 X 107-7 In Tayg )N*”°? (Th — Te) + 


Ce = Tr) +4 Ged, _ TN, (4) 
Table 2: MLI properties in two blankets [21]. raa = TShroud~Tiut outer (5) 

ML inner ML outer Exterior shroud &€MLI outer 
MLI Layer Density, N, (cm ') 17.7 17.7 The degradation factor, D,, is adjusted 
MLI Layers, N 17 wi to match the overall heat flow, 48.2 W. 
Emissivity, € 0.04 0.04 0.05 Results show (Figure 9) that total external 
Degradation factor, D, 2.7216 wall heat flow is nearly constant 
Interstitial Gas Pressure, P. (torr) 1.x10° throughout the simulations. The heat flux 
Shroud Temperature, K 350 values for this tank wall boundary 
Modified Lockheed | C,=2.4x107; C,=4.944x107°: condition vary less than 5% over the entire 
Equation constants C,=1.46x10* surface, while wall temperature varies by 

~25 K. 


E. Initial Conditions for K-Site Tank Simulations 

The original experimental report [1] noted that, “the starting conditions had a significant impact on the transient 
initial pressure rise rate.” Two starting conditions, ‘Steady Boil-Off’ and ‘Isothermal’, were tested at the 29% and 
83% fill levels, but at the 49% fill level only Steady Boil-Off initial conditions were tested. Figure 5 shows the 
experimentally observed pressure difference between the two initial conditions for 83% fill—within the first hour, 
an ~8 kPa pressure difference develops between the two conditions; the pressure difference for 29% fill is similar 
[1]. 

How were these different starting conditions achieved experimentally? Paraphrasing Ref. [1], at the conclusion 
of the previous test, the tank was slowly drained to just above the desired fill level, and tank pressure was reduced to 
atmospheric pressure (103 kPa). This venting “induces substantial bulk boiling of the fluid that initially produces 
nearly isothermal conditions in the tank.” The Isothermal initial condition started from this condition, in particular, 

“shortly after tank pressure is reduced to 

0.4 atmospheric pressure by venting. As soon as the 
tank lid temperature reaches its minimum value 
(~23 K), the test is begun by closing the vent line 
valves.” For the Steady Boil-Off starting condition, 


E = tank venting was maintained for 4 hours or more 
g until the liquid surface-to-tank lid temperature 
00 +-----Sgse————— wee +--+ ------- gradient and boil-off rate stabilized. 
c eae . . . 
Bis Boiling is a high heat transfer mechanism— 
£ me anerace (NIST Date) much higher than thermal conduction. Bulk boiling 
g O Interface (Reynolds Data) 2 : < . 
ge 02 > Osarientlbiaae with vapor bubbles moving through the liquid 
2 os —sothermal Case Initial T Profile would tend to create uniform liquid temperatures; 
8 —Lower T_Bulk Initial T Profile vapor bubbles boiling at the wall and moving to the 
§ -o4 sw Steady BolhOfInidal:T Rrofile surface through colder fluid (below T,,,(p)) would 
> — -Interface : 5 F 
be contract as condensation heated the fluid while 


those moving through warmer fluid (above T,,,(p)) 
would expand as evaporation cooled the fluid. 

What deviations from truly isothermal 
conditions might be expected in the Isothermal 
starting condition? With a pressure drop, 
evaporation will depress the fluid temperature near 


20.0 20.2 20.4 20.6 20.8 21.0 21.2 21.4 21.6 21.8 22.0 
Temperature (K) 
Figure 7: Initial temperature profiles for simulations. 
The initial pressure transient is very sensitive to the 
initial temperature profile near the interface. 
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Table 3: Simulations, their grids, and initialization. Initial velocity zero. Experimentally measured temperature 
profile used for initialization in vapor phase and near the interface; bulk liquid temperature elsewhere. 


Time Interface Wall Normal Initial Bulk 
Case Name ID Grid Cells | Step Normal Spacing Liquid 
(ms) | Spacing (nm) (mm) Temperature (K) 
Tsothermal-Fine K8 4 x Finer 61,400 1 0.5 0.25 20.3 
Isothermal K10 Baseline 14,800 10 0.5 0.5 20.3 
Steady Boil-Off | K12 | Improved Baseline | 15,100 10 0.5 0.5 20.3 
Lower Tau 06 2 x Coarser 9,800 20 2.0 0.5 20.25 
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Figure 8: Isotherms help reveal heat flows in the K-Site tank 
simulation. Heat flows shown by orange arrows and fluid motion shown 
by white arrows. Heat pathways indicated by letters referenced to Figures 
9 and 10. Heat flows change over time and reach equilibrium after O(10) 


hours. 
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the interface as heat is consumed for 


evaporation. Typically, boiling 
requires a degree or two. of 
superheating above saturation 


temperature, so boiling will end before 
true isothermal conditions are reached, 
and liquid near the walls will be above 
Tya(p). A smaller effect is pressure 
head. For the 49% fill LH2 case, the 
pressure at the tank bottom is 650 Pa 


higher than the interface, 
corresponding to a saturation 
temperature 0.022 K — warmer. 


Certainly, isothermal conditions only 
hold approximately. Further, from 
available experimental data [17], there 
is only data for one temperature sensor 
in the bulk liquid. Its initial value is 
20.256 K and it remains below 20.354 
K for over an hour. 

Numerical simulations indicate that 
the initial pressure transient (Figure 5) 
is very sensitive to the initial 
temperature profile near the interface 
(Figure 7). The Isothermal case was 
initialized with the initial temperature 
profile from experimental data [17] in 
the vapor phase and near the interface, 
plus 7);, = 20.3 K (Figure 7) in the bulk 
liquid. The Lower Tg, case changes 
only Tig = 20.25 K. To simulate 
starting conditions, the Steady Boil- 
Off simulation was run with several 
hours of venting at p = 103.632 kPa 
before self-pressurization, 80 Tinerface = 
Tsa(p=103.632 kPa) = 20.354 K. But 
note the complex profile in Figure 7 
with temperature depressed at the 
interface. 


F. Simulation Operational Details: 
Starting Simulations and 
Monitoring Time Step Convergence 
Baseline and fine grid simulations 
are hard to start even from good initial 
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conditions; coarse grids are more forgiving. Initial iterations at small time steps will start these grids: 300 iterations 
at 10s, 300 at 10s, 3000 at 10s, 3000 at 2x10™s, 3000 at 4x10°s. 

The interface temperature must be close to saturation conditions, Tya/(Dinitiqi), Or the simulation will start with 
large amounts of evaporation or condensation as the solution moves the interface to saturation conditions; solutions 


may not tolerate this shock. 


To detect potential problems, solutions are monitored for global mass conservation, and global total energy 
conservation. It is easy to miss convergence at a time step since it is buried in Fluent’s output. An awk script scans 
output to detect sporadic, malignant convergence failure in individual time steps. During the simulation, longer time 
steps are tested to see if they are possible while maintaining convergence. Table 3 gives the highest time steps used. 


IV. Results: Observed Phenomena and Behavior 


The original self-pressurization experiment report noted that “the mode of heat transfer is complex and is the 
greatest factor controlling the pressure rise rate.” [1, p. 9]. The simulations have been analyzed to identify and 
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Figure 9: Global heat flows within the tank simulation 
with time. Heat flows are identified by letter, 
corresponding to flows in Figure 8. Isothermal case. 
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Figure 10: Heat flow tangential to the wall as measured 
by various pairs of temperature sensors. The aluminum 
walls are very conductive: the ratio of thermal 
conductivities for GH2:LH2:Al is 1:6:1180. Letters 
correspond to Figure 8. Isothermal case. 
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quantify heat transfer mechanisms as well as fluid 
behavior. 


A. Heat Transfer Paths 

Figure 8 shows the heat transfer paths identified 
in the simulations. Important paths include heat 
transfer through the MLI (denoted A), substantial 
heat conduction tangentially through the thin—yet 
highly conductive—tank walls (denoted W,, W2, 
W3), much more heat transfer from the wall into the 
liquid than the vapor (B, C), natural convection 
along the tank walls, and condensation (D) and 
evaporation (E) at the liquid/vapor interface. 
Smaller heat flows include conduction (F) through 
the vapor phase towards the interface, and 
convection, turbulent conduction (G) through the 
liquid phase towards the interface. 

Figures 9 and 10 show measurements of these 
heat flows over time marked with letters. 
Importantly, these heat flow rates change with time, 
and only reach equilibrium after O(10) hours. 

There is so much heat flow through the walls 
into the liquid, as noted above, because aluminum 
has relatively high thermal conductivity. Also, the 
LH2 can absorb much more heat than GH2 for the 
same temperature rise. 

I. Measuring Heat Flows 

How do we quantify and calculate these heat 
flows? The Fluent solution includes boundary heat 
fluxes that are integrated over the wall inner and 
outer surfaces (heat flows A, B, C) every 100 time 
steps. For condensation and evaporation (denoted 
D, E), mass transfer is integrated over the volume. 

Temperatures at specific locations are logged 
through the simulation and used to estimate the 
conductive heat flow, Eq. (6), and how it changes 
through the simulation. 


Q=k 2a (6) 
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Figure 11: Temperature contours (AT=0.05 K) 
contact point and (right) velocity vectors in the surrounding region. 
Red cross circles indicate temperature monitor locations for calculating 
heat flux and estimating heat flow. The temperature monitor pairs W, 
and W; are used to calculate heat flow through the wall near the contact 


point. The F temperature pair calculates vapor conductive heat flow. 
m/s. 
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Figure 12: Heat Flows at the interface. 
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For example, heat flow tangential to 
the wall near the contact point (denoted 
W,, W2) is calculated from the 
temperature monitor pairs shown in 
Figure 11—plus an estimate of the 
wall cross-sectional area, A. Note that 
a point estimate of heat flux is 
extrapolated to a heat flow for an area, 


ee Layer 


Heat conduction in the vapor 
phase is measured with Eq. (6) and 
many pairs of temperature monitors, 
like those shown in Figure 11. From 
many measurements, the typical vapor 
phase heat flow is 1.5-3.5 W at a fluid 
level—small for such a large fluid 
level area. A separate simulation 
(Steady Boil-Off) on a_ similar 
size/resolution grid gives similar heat 
flux values. 


(left) near the 


B. Pressure Evolution, Interface 
Temperature, and Heat Fluxes into 
the Interface 

Vapor pressure is not only the 
saturation pressure, Psa(T interface), at the 
liquid/vapor interface, but it reflects an average 
temperature of the interface, Tinerjace. Deviations from a 
single value of Tinterface—cold or hot, liquid or gas, 
impinging on the interface—results in evaporation in 
warm areas and condensation in cold areas that pushes 
the interface temperature towards a single interface 
value. 

Further, heat fluxes into the interface supply or 
remove the heat energy required for evaporation or 
condensation—influencing pressure evolution. Figure 12 


demonstrates interfacial heat flows 
conceptually, and Figure 13 gives 
interfacial heat flow estimates from a rake 
of five temperature monitors spanning 1 
cm of the interface, using Eq. (6) and 
extrapolated to the entire interface. After 
an initial transient, the vapor to interface 
conductive heat flow reaches a value of 5 
W into the interface. Convective heat 
flow is hard to measure, and could be 
small. Many temperature monitor pairs in 
the gas suggest similar gas phase heat 
fluxes. 

Similarly, the interface to liquid 
conductive heat flow is small. Convective 


"Gas at Interface 
“wr"Bulk Gas Probes 2-8 
"Bulk Gas Probes 1-2 


Liquid at Interface 


18 


Most 


heat transfer is more likely here. 
evaporation occurs near where 
meniscus meets the wall. 


A 


Figure 13: Estimated heat flows at and near the interface. 
rake of five temperature monitors in the simulation provides 
estimates of heat flow using Eq. (6). Note that a point value of heat 


the 


flux is extrapolated to the entire interface. Heat flow is vertically 


downwards. Bulk gas monitors 1, 2, 8 suggest similar values. 
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Figure 14: Initial pressure evolution in experiments 
and simulations. The Lower Tp, case at 49% fill 
matches experimental data best; it uses an 0.05 K lower 
bulk liquid temperature, which, presumably, reduces 
heat conduction from liquid to the interface, and slows 
interface temperature increase. The closeness of the 


match is much less relevant than the demonstrated 
Experimental 


temperature _ sensitivity. 
accuracy +0.01 kPa. 
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Figure 15: Pressure evolution during experiment 
and simulations. The Isothermal case at 49% fill 
displays the pressure evolution that would be expected 
for Isothermal initial conditions, but no experimental 
data exists for comparison. The Lower Tau, case 
closely matches the Steady Boil-Off experimental 
conditions at 49% fill level until 8 hours when it under 
predicts pressure. The Lower Ty, case uses a 0.05 K 
lower bulk liquid temperature (Figure 7). 


C. Natural Convection Up the Walls 

Figure 8 shows that buoyantly driven flow (natural 
convection) occurs in the liquid along the wall, and to a 
lesser extent in the vapor near the wall. Heat transfer from 
the wall to the liquid creates a thin thermal boundary layer 
of warm liquid that is buoyantly driven up the wall. This 
velocity creates a thin momentum boundary layer (Figure 
11). The buoyantly driven liquid accumulates in a stable 
layer of warm liquid near the interface. Eventually, this 
layer is warm and thick enough to buoyantly decelerate the 
upwelling flow along the wall, and the flow cannot 
penetrate the layer. Figure 11 demonstrates this buoyancy 
barrier in contrast to Figure 18 where the barrier has not 
formed early in the simulation. This phenomena 
contributes to the complex heat transfer process. 


D. Less Relevant Phenomena and Behavior 

I, Internal Gravity Waves 

Obviously, the interface can have waves if it is 
perturbed, but internal gravity waves would also occur 
within the density stratified vapor and liquid phases. 
These waves are observed numerically at very short time 
steps. Internal gravity waves should not significantly 
influence mixing unless the waves break. Breaking near 
the rim of the manhole access is likely. Similarly, wave 
breaking around instrumentation rakes is likely. 


V. Results and Comparison with Experimental 
Data 


A. Initial Pressure Evolution 

As noted in Section HI-E, the initial pressure evolution 
is very sensitive to initial conditions: Isothermal and 
Steady Boil-Off. Figure 5 shows a pressure difference of 
~8 kPa develops over the first two hours in the 83% fill 
case for the two different initial conditions. The 29% fill 
case is similar. 

Several initial temperature profiles (Figure 7) were 
tested in simulations for the 49% fill level, and Figure 14 
shows three results. The Steady Boil-Off case tried to 
initialize with several hours of evaporative venting, but 
with only a modest improvement. Yet, the best agreement 
with experimental data is in the Lower Tay, case (Table 
3), where a 0.05 K lower bulk liquid temperature was 
used—20.25 K instead of 20.3 K. From available 
experimental data [17], there is only one temperature 
sensor in the bulk liquid, with an initial value of 20.256 K. 
Of course, all these temperature differences are well 
within the experimental temperature sensor accuracy 
range. 

How do we interpret this behavior? The initial 
pressure rise, P.at(Tinterface), reflects the temperature rise of 
the interface and the local heat flows that determine it. A 
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Figure 16: Comparison of temperature evolution at five experimental and 
simulated fluid sensors. The Steady Boil-Off case shows good agreement except 
at the top sensor in the lid. The tank lid is simulated as a solid piece neglecting 
the reduced thermal contact conductance at the bolted joint. The simulation 
temperatures for negative time are solving venting conditions for the initial 
conditions. Experimental temperature sensors calibrated to +0.1 K. 
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Figure 17: Comparison of temperature evolution at six experimental and 
simulated wall sensors. The Steady Boil-Off case shows good agreement except 
for sensors near the tank top and lid. The tank lid is simulated as a solid piece 
neglecting the reduced thermal contact conductance at the bolted joint. The 
simulation temperatures for negative time are solving venting conditions for the 
initial conditions. Experimental temperature sensors calibrated to +0.6 K 
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lower bulk fluid temperature 
reduces the heat flow to the 
interface and slows the 
interface temperature rise. 

In conclusion, the closeness 
of the match between 
experimental and Lower Tpux 
case data is much less relevant 
than the demonstrated 
temperature sensitivity. 


B. Pressure Evolution 

Figure 15 compares the 
pressure evolution during the 
experiment with two 
computational results. The 
available experimental data [1] 
[17] at 49% fill is for Steady 
Boil-Off initial conditions. 
The Lower Tpux case agrees 
with experimental results up to 
8 hours when it under predicts 
the pressure rise. This case 
uses a 0.05 K lower bulk liquid 
temperature. 

The Isothermal case data 
display the behavior that would 
be expected for Isothermal 
initial conditions, but no 
experimental data exists for 
comparison. 


C. Fluid and Wall 
Temperature Evolution 

The agreement between 
experimental and simulation 
fluid and wall temperature 
sensors is good except near the 
tank lid and top as shown in 
Figure 16 and Figure 17. This 
difference is attributed to 
treating the man access lid 
(Figure 2) as a solid piece of 
aluminum and neglecting the 
reduced thermal contact 
conductance of the bolted 
joint. Studies of heat transfer 
through bolted joints indicate 
reduced heat transfer [18]. 
This issue should be corrected 
in a future simulation of the K- 
Site tank. 


VI. Numerical Observations 


In performing these simulations, several numerical subtleties were observed that deserve mention as they 
improved these simulations and are important for future simulations. 


A. Mass Loss in the Simulations and Its Resolution 
Early K-Site tank simulations showed a slow fluid 
mass loss of ~500 grams over several hours, 
corresponding to an interface level drop of several 
millimeters. This mass loss was traced to the mass 


Unphysical transfer UDF, Fluent’s Define Mass Transfer(), and 
Condensation in 7 c. 


which grid cells included mass transfer, m. If mass 
transfer was allowed where 0.05 < c < 0.95—between 
the red interface boundaries in Figure 18—then mass 
was conserved 10° times better than in the traditional 
region—between the black boundaries. The 
interpretation is that an unphysical condition exists 
when condensation is specified in a cell containing all 
liquid and no vapor. Further, Fluent’s mass transfer 
implementation assumes physical consistency, and the 
4 : . Theory Manual [9, p. 608] shows one way to filter out 


1| fs « : 7 - these inconsistencies. Eq. (3) excludes condensation in 
ths id ~ o liquid filled cells and evaporation in vapor filled cells. 
ry a ee = é re 
Heat s * - oe a - ~ 
Flow 7 et = ~ s r au 3 
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Figure 18: Interface region near the contact point Figure 19: Spurious velocity vectors along the 
showing velocity vectors (below) with inset for  liquid/vapor interface (red lines)—'Where does this 
(middle) P,,.-P plus thresholds (original in black; momentum come from?' Surface tension modeling 
0.05<c<0.95 in red) for calculating mass transfer with with VOF results in these spurious velocities. From 
inset for detail (top). | Unphysically specifying K5 simulation at 17820 s; velocity range [0, 0.157] 
condensation in liquid filled cells resulted in mass loss. m/s. 


B. Surface Tension Modeling 

Initial simulations found spurious velocities along the liquid/vapor interface as shown in Figure 19. The 
operative question was: Where is this momentum coming from? This issue was traced to surface tension modeling, 
which the literature [19] indicates is a known issue with VOF. Surface tension modeling was not included in 
subsequent simulations, and the behavior of the interface was much more plausible as shown in Figure 11 and 
Figure 18. Since including surface tension would only resolve the meniscus (~2 mm high) at the contact point— 
which is poorly resolved by the grids anyway—this is an appropriate modeling choice here. 


C. Evaporation at the Meniscus/Contact Point 
Another phenomena contributing to the complexity of this simulation is evaporation at the meniscus where the 
liquid/vapor interface meets the wall. As measured in Figure 10 and demonstrated in Figure 18, there is a 


12 
American Institute of Aeronautics and Astronautics 


substantial heat flow (10-15 W) through the aluminum wall at the contact point. Further, simulations show a point 
of strong evaporation on the wall at the contact point as shown in Figure 18 and Figure 20. 

vee can the heat in the wall go? Above the contact point, heat transfer is limited by conduction into a very 
non-conductive gas—GH2 is 1180 times less 
conductive than aluminum. Similarly, below the 
contact point, heat transfer is limited by heat 
conduction into a non-conductive liquid—LH2 is 200 
times less conductive than aluminum. Yet 
evaporative heat transfer is possible in between at the 
thinnest part of the meniscus: the wall heats the 
liquid on one side of the film, and it evaporates on the 
other. The liquid’s thermal conductivity doesn’t slow 
heat transfer when the meniscus is so_ thin. 
Independent models indicate substantial heat flow 
Simulation and evaporation in the thinnest art of the meniscus. 
ae: Some experiments have been done with a sub-grid 
model for this additional evaporation source. 

Solutions to the Young-Laplace equations [20] 
indicate a meniscus height of ~2 mm for contact 
angles less than 10 degrees—3-4 grid cells on the 
baseline Isothermal grid as shown in Figure 20. 
Despite the highest resolution of the grid, the 
simulation may not resolve this evaporation. 


Figure 20: Comparison of VOF interface and 
actual meniscus shape. The meniscus height is ~2 
mm for contact angles less than 10 degrees. There is a 
point of high evaporation in the simulation, and high 
evaporation is expected at the thinnest point in the 
meniscus. Both points are enabled by high heat flux 
through the wall. 


VIL. Conclusion 


This paper documents the successful implementation of detailed tank geometry and MLI insulation into a 
cryogenic hydrogen tank simulation. Self-pressurization has been carefully studied and simulated. In particular, 
heat flows within the tank have been identified and quantified based on multiple simulations. Further, the sensitivity 
of the initial pressure rise to the initial temperature profile near the interface has been identified. Numerical issues 
have been identified and remedied, and they will improve future cryogenic fluid simulations. These issues include 
surface tension modeling, improved representation of saturation pressure, and identification and repair of mass loss 
in the simulation. Clearly, both experimental measurement accuracy and uncertainty in saturation pressure, P,,(T), 
complicate simulation of pressure evolution. Better calibration of experimental temperature sensors would be 
helpful, particularly at the interface and in the LH2. That velocity variations are introduced at the interface is a 
concer for the VOF method. Real gas effects should be accounted for in future simulations. Thermal contact 
conductance at bolted joints should also be included. 
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